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^ Abstract 
> 

• ^ We propose a novel algorithm, called REGGAE, for the generation of momenta 

of a given sample of particle masses, evenly distributed in Lorentz invariant phase 
^ space and obeying energy and momentum conservation. In comparison to other 

existing algorithms, REGGAE is designed for the use in multiparticle production in 
hadronic and nuclear collisions where many hadrons are produced and a large part 
of the available energy is stored in the form of their masses. The algorithm uses a 
loop simulating multiple collisions which lead to production of configurations with 
reasonably large weights. 
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Program summciry 



Program title: REGGAE (REscatterig-after-Genbod GenerAtor of Events) 
Catalogue identifier: 

Program summary URL: http://www.fpv.umb.sk/~tomasik/reggae 

Program obtainable from: http://www.fpv.umb.sk/~tomasik/rcggae 

RAM required to execute with typical data: This depends on the number of 

particles which are generated. For 10 particles like in the attached example it 

requires about 120 kB. 

Number of processors used: 1 

Computer (s) for which the program has been designed: PC Pentium 4, though 
no particular tuning for this machine was performed. 

Operating system(s) for which the program has been designed: Originally de- 
signed on Linux PG with g-|--|-, but it has been compiled and ran successfully 
on OS X and MS Windows with Microsoft Visual C++ 2008 Express Edition, 
as well. 

Programming language: G++ 
Size of the package: 12 KB 
Distribution format: zipped archive 

Number of lines in distributed program, including test data etc.: 1468 
Number of bytes in distributed program, including test data etc.: 52 KB 
Nature of physical problem: The task is to generate momenta of a sample of 
particles with given masses which obey energy and momentum conservation. 
Generated samples should be evenly distributed in the available Lorentz in- 
variant phase space. 

Method of solving the problem: In general, the algorithm works in two steps. 
First, all momenta are generated with the GENBOD algorithm. There, particle 
production is modelled as a sequence of two-body decays of heavy resonances. 
After all momenta are generated this way, they are reshuffled. Each particle 
undergoes a collision with some other partner such that in the pair centre 
of mass system the new directions of momenta are distributed isotropically. 
After each particle collides only a few times, the momenta are distributed 
evenly across the whole available phase space. Starting with GENBOD is not 
essential for the procedure but it improves the performance. 
Typical running time: This depends on the number of particles and number 
of events one wants to generate. On a LINUX PC with 2 GHz processor, 
generation of 1000 events with 10 particles each takes about 3 s. 



* Supported in parts by ITMS:26220120007 (Slovakia), MSM 6840770039, and 
LC 07048 (Czech Republic). 
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1 Introduction 



A frequent task encountered in simulations of particle production is to generate 
a sample of many particles which exactly conserve energy and momentum. 
The events generated in the procedure should be uniformly distributed in the 
Lorentz invariant phase space (LIPS). This is a standard textbook problem 
in the case of two or three-body emission, but becomes increasingly involved 
as the number of particles grows. Such situations may also become relevant 
now with the high multiplicity data from proton-proton as well as Pb-|-Pb 
collisions at the LHC Algorithms have been developed based on the formulas 
for integration of the LIPS which can be used for such a task. Among them 
we find GENBOD based on the treatment by sequential decay [1,2]. 

Originally, the motivation for the development of these algorithms was in the 
need to simulate the Lorentz invariant phase space for Monte Carlo integra- 
tion of complicated matrix elements. The problem of simple algorithms like 
GENBOD is that for a large number of particles it generates too often con- 
figurations which are very unlikely and have tiny weight of their contribution 
to the phase space integral. If the algorithm is used for event generation, 
such configurations are accepted with very low probability and this makes the 
simulation rather inefficient. 

Improved algorithms are available on the market like RAMBO [3] and NUP- 
HAZ [4]. Those, however, are mostly aimed for the use in high energy pro- 
cesses with rather large amount of energy in the form of kinetic energy, i.e. 
the masses small in comparison to the energy. The effectiveness drops fast if 
the amount of energy contained in the masses grows considerably. This is a 
serious drawback which makes these algorithms unusable for the simulation 
of multiparticle production in nuclear collisions where masses are of the order 
of, or larger than, the momenta. Also, high multiplicity proton-proton data 
are beyond reasonable applicability of these tools. 

When working in this regime a different kind of approach is sometimes used. 
One rather starts with generating the particles from thermal distribution. The 
momentum of the last particle is then calculated from momentum conservation 
and the energies of all particles are rescaled so as to match the total energy [5] . 
Another possibility is to calculate the energies and momenta of the last pair 
of particles [6] . These algorithms are not guaranteed to always work and one 
may have to repeat them in order to generate a usable configuration. Moreover, 
the particles with the calculated momenta may lie far out off the region where 
other thermal ones are concentrated. The chance of the appearance of such 
problems grows with the number of particles. 

Recently, a new generator was reported in [7] but no details are available to 
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us about its construction and performance. 

In our approach we first generate n-body events which do conserve energy 
and momentum but are not guaranteed to fill the Lorentz invariant phase 
space uniformly. Then, the generated momenta are reshuffled in (virtual) two- 
body collisions. A few collisions per particle lead the system to the most 
likely configuration. In this way we solve the problem of generating unlikely 
configurations: our algorithm samples the space of all possible configurations 
according to their probability of appearance. This leads to uniform distribution 
of events within the available LIPS. 

Details of the algorithm are explained in Section 2. This is followed in Sec- 
tion 3 by a demonstration of the approach of events generated by REGGAE 
to equilibrium after increasing the number of collisions. To this end we esti- 
mate the information entropy of the generated configuration. In Section 4 we 
compare the results of REGGAE with RAMBO, NUPHAZ, and the original 
GENBOD algorithm. A short manual for the C-|— I- routines which accompany 
the paper is given in Section 5. 



2 REGGAE: the algorithm 



We must first generate hadron momenta so that they will conserve the fixed 
total momentum and energy. Any procedure fulfilling this simple requirement 
can be used, because the final configuration of momenta is formed later in the 
rescattering part. A good choice of the initial generator can, however, make the 
algorithm more efficient, as less improvement is required from the subsequent 
rescattering part. Our procedure starts with GENBOD [2]. 

We want to distribute energy E and momentum P among n particles with 
masses mi, m2, . . .m„. We can boost into the frame where P vanishes and 
the energy assumes the value E*. After the generation is complete we boost 
back to the original frame. The trick is to formally treat the multiparticle 
generation as a sequence of two-body resonance decays. First, an auxilliary 
resonance with mass E* decays into a particle with mass m„ and a resonance 
with mass M^-i- Then M„_i goes into m„_i and M„_2 etc. Directions of the 
momenta of decay products are random up to momentum conservation. The 
essence of the procedure is in determining the masses of auxiliary resonances. 
They are chosen randomly but they must fulfill the inequalities 

rui + Mi_i <Mi< Mi+i - rrii+i , (1) 

where Mi = mi and M„ = E* . The condition is fulfilled if we generate Mj's 
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via 

Mi^j2^k + xjE*-J2mi] , (2) 

k=l \ 1=1 ) 

where Xi is random variable from the interval (0, 1). Condition (1) is satisfied 
if we choose Xj's so that 

< < X4 < • • • < x^-x ■ (3) 
Details of the procedure can be found in [2] . 

As wc explained in the introduction, momentum configurations generated by 
GENBOD often come out with very small weight. In practice, this means that 
such an event is extremely unlikely to appear in Nature. Therefore, in the next 
step we reshuffle the momenta in order to achieve a more likely configuration. 

In the routine, the momenta of all particles are stored in an array. We randomly 
choose a distance d between positions within the array. Then each particle 
collides with a partner that is stored by d positions further down the array. 
Particles in the last d positions collide with partners in the beginning of the 
array. In every turn each particle collides twice. 

The coUisions are simple s-wave scatterings. We always Lorentz boost into 
the center of mass system of the colliding pair, generate new directions of 
the momenta with isotropic distribution, and then boost back to the original 
frame. For a small number of particles 6 collisions per particle is enough to 
reach a distribution with large weight. This number becomes larger for large 
numbers of particles. 



3 The approach to most likely configurations 



For very large number of rescatterings our procedure leads to uniform filling 
of the available Lorentz invariant phase space. The first step to prove this is to 
construct the distribution of outgoing particles in two-particle collisions. The 
distribution is defined as the probability dP of finding the outgoing momenta 
(Pi,P2) in a given elementary volume d^pid^p2 of the two-particle momentum 
space, given the incoming momenta (51,52)- the distribution is isotropic in 
the center of mass system, we have in a laboratory system 



where and are the total four-momenta of incoming and outgoing par- 
ticles, respectively, (£^i,£^2) are the energies of outgoing particles, and 
and /i^ is the sum and the difference of the masses of particles squared. 
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Af ^ = mf + m\ and [x^ = ml — m|. The crucial observation is that the distri- 
bution in (4) normahzes not only in outgoing momenta but also in incoming 
ones, i.e. it satisfies 

where (61,62) are the energies of incoming particles. This is seen immediately 
after replacing in the expression in front of (5-f unction by p^, since (5) 
then reduces to the normalization in outgoing momenta with (Pi,P2) renamed 
to {qi,q2) s-nd vice versa. Thanks to (5), collisions with the distribution of 
outgoing particles (4) produce an equilibrium distribution of the system of n 
particles that is uniform in LIPS 

dp^Cp5\p,,,-P)ll^, (6) 

1=1 * 

where pf^t is the total four-momentum of the system, P^^ is the value as- 
signed to Ptot and Cp is a normalization constant depending only on P'^. 
Indeed, if we view dP for an arbitrary pair of particles as a two-particle block 
of a block-diagonal n-particle transfer matrix, it follows from (5) that (6) 
is an eigenvector of the transfer matrix corresponding to the eigenvalue 1. 
Then, Perron-Probenius theorem guarantees that it is the only eigenvector 
with that property. (The last claim holds if the particle momenta can change 
from any values to any other values consistent with the conservation of energy- 
momentum in a finite number of collisions, which seems plausible enough.) The 
state of the system described by the eigenvector of the transfer matrix with 
the eigenvalue 1 is obviously an equilibrium state, since it does not change 
in the evolution of the system. Another consequence of the Perron-Frobenius 
theorem is that the system evolves towards it from any nonequilibrium initial 
state. 

While this proves that our procedure will eventually lead to the uniform filling 
of the available LIPS, we also want to demonstrate how fast this happens, i.e., 
how fast the system approaches the most likely (equilibrium) configurations 
when starting from unhkely ones generated by GENBOD. In algorithms based 
on Monte Carlo integration of the phase space usually a weight is determined 
for each generated configuration. This weight could be used as a measure of 
likelihood of a given configuration, however, it depends on integration variables 
and it is not clear that we can define it uniquely in our case. 

We thus use a different approach, with information entropy as a measure of 
likelihood, defined not for a single configuration, but for the whole set of 
Ne configurations (events), each with n particles. We assume that the mo- 
menta of all nNg particles are distributed according to the same underlying 
single-particle probability density distribution p{p). The information entropy 
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is defined as 

S^- l^pip) In pip) d^p (7) 
where the integration runs through all accessible momentum space E. 

We want to show that the rescattering part of the algorithm changes this en- 
tropy in such a way that it will grow from some value corresponding to unlikely 
configurations to higher values for more likely configurations and eventually it 
will saturate at equilibrium. As can be seen, we formulate this calculation for 
a distribution on a three-dimensional momentum space. The distribution p(p) 
is generally unknown. We estimate it from the sample of momenta from the 
set of Ne configurations: first the available three-momentum space is divided 
into N elementary cells of equal volume AV^, then the probabihty density p(p) 
in a cell AVi is estimated by pi which is given as 

where is the number of particles which fall in the cell AVi. 
Then the entropy is estimated as 

N 

S^-Y.iPi^^Pi)^Vi^ (9) 
1=1 

where is large and AVi ^ 1/-^ small. 

At equilibrium the information entropy can be calculated theoretically in the 
limit of large n, since p{p) is known in this case. It is related to the single- 
particle energy distribution Pe{E) (called here the LIPS-Boltzmann distri- 
bution), which is derived from the uniform distribution in LIPS using the 
Darwin-Fowler method and is given by [8] 

PE(E)dE^AfVE^ - m^exp (--^) dE, (10) 

where E = y/p'^ + is the total energy of a single particle, T is the LIPS- 
temperature, m is the mass of the particle and ks the Boltzmann's constant. 

is a normalization constant, and Ki (m/ZcsT) is the modified Bessel function. 
Note that Pe{E) differs from the canonical Boltzmann distribution in two im- 
portant aspects: the canonical distribution has an extra E in front of the expo- 
nential and the canonical temperature is different from the LIPS-temperature 
T. The differences come essentially from the fact that unlike the canonical rel- 
ativistic Boltzmann gas which thermalizes via collisions both in configuration 



7 



and momentum space, our system thermalizes via collisions which take place 
only in the momentum space. The LIPS-temperature T is implicitly defined 
by the equation for the mean value of the single-particle energy 



(E) = / EpE{E)dE 



— m 



K^jm/kBT) ^ E*_ 
Ki(m/A;sT) n 



(11) 



Then, the equilibrium probability density p(p) is found from the LIPS-Boltzmann 
distribution 



The equilibrium information entropy is then found by substituting p{p) into 
the definition of S in Eq. (7) and integrating. For configurations of n = 100 
pions (m = 139 MeV) of the total energy E* = 50 GeV we get the LIPS- 
temperature ksT — 0.20688 GeV and the equilibrium information entropy 
prediction S = 0.6724. Our expectation is that the entropy estimated from 
Eq. (9) will grow with the number of collisions Nc and for some A^^^ it will 
saturate at the equilibrium value which should coincide with the theoretical 
prediction. We show the result of this calculation in Figure 1. 

The entropy saturates at the predicted equilibrium value after 12 colisions in 
this particular case. For a smaller n, we found that the entropy saturates after 
a smaller number of collisions. Our findings are in line with studies of kinetic 
theory. For example, in [9] it has been shown that a few collisions are enough 
to drive the system towards equilibrium distribution. 

This indicates that our algorithm can generate relevant momentum configu- 
rations not only in principle, but it can do so in a reasonably small number of 
steps. Further support that this is indeed so will come from the next section 
where we apply REGGAE to Monte Carlo integration. 

Note here also, that the rescattering algorithm would arrange the momenta 
correctly even if we did not start with GENBOD. We have checked this by 
starting with 15 equal momenta all in the same direction and another 15 ones 
of equal size and opposite direction. Rescattering succeeded in arranging the 
momenta but needed more steps to reach saturation. This shows that starting 
with GENBOD is not essential, but improves the effectivness. Calculations 
using GENBOD were indeed faster. 



pip) 




(12) 
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Fig. 1. Evolution of the information entropy estimates with the number of colhsions 
Nc for 10^ REGGAE-generated configurations of n = 100 pions (m = 139 MeV) 
of the total energy E* = 50 GeV. The red line shows the theoretical equilibrium 
value S = 0.6724, the blue curve shows the estimates of Eq. 9 based on REGGAE 
configurations. The statistical errors of the estimates as < 0.0001. 

4 Comparison with other algorithms 



One might argue that the information entropy saturation at the equilibrium 
value demonstrated in the previous section only shows that the single-particle 
momentum distribution has LIPS-thermalized while it is not clear that many- 
particle momentum distributions have done the same. To provide further sup- 
port for a fast, few-step LIPS-thermalization, we applied our algorithm to 
Monte Carlo integration over the phase space, a procedure sensitive to many- 
particle distributions. In this section we also compare REGGAE with other 
routines available on the market. 



Generally, the infinitesimal element of the LIPS for n particles is 

-Pn) (13) 

The general task is to calculate integral over available phase space 



= (2vr)32i^;, (2vr)32E, ' ' ' Jt^Y^J i^-Pi-P^-Ps' 



f{{p^})d<!>r. 



{Pi} = {Pl,P2,P3, ■■■,Pn) 



(14) 
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where f{{pi}) is some function of the momenta. Algorithms for Monte Carlo 
integration generate configurations of momenta, calculate their weights Wj and 
then approximate the integral with (sample mean method) 



' Ef=i^/(fa}. 



$n, (15) 



where Wj = Wj/wmax is the ratio of the weight of the configuration to the max- 
imum weight calculated or found empirically in the large set of configurations 
and the sum runs over N different momentum configurations {pi}- 

If the configurations correspond to real events, they have w', = 1 and 



1 ^ 

T^E/(M.)| (16) 



We have tested the sum of Eq. (16) calculated with events generated by our 
algorithm against calculations by RAMBO [3], NUPHAZ [4], and GENBOD 
[2] . In case of GENBOD we accepted the generated events with the probabil- 
ity determined by the ratio of the event weight and the maximum theoretical 
event weight. This we call weighted GENBOD (wGENBOD). It should be dis- 
tinguished from the unweighted GENBOD (uGENBOD) where we ignore the 
weights, which is incorrect and we show it here for the sake of demonstration. 
We have performed the tests on a variety of functions. Here we present just 
one example. 

We integrated a function of five four-momenta 

Mpi,P2,p3,P4,P5) = m^+pIpI ^ ^ 

with — 25 GeV^. The phase space was given by the total four- momentum 
P = (100, 0, 0, 0) GeV which was distributed among particles with the mass 
of m = 1 GeV. 

In Table 1 we compare the results of the five generators together with the 
times of computation of 10^ events on AMD Athlon(tm) 64 X2 Dual Core 
Processor (speed: 3 GHz, RAM: 2 GB, OS: Ubuntu Linux, kernel: 2.6.31-22). 
We show the mean values of /s, i.e. only the part of formula (16) within 
the parentheses. The algorithms, including REGGAE, yield consistent results 
except for the unweighted GENBOD, which was expected. Although the result 
of REGGAE is fine, its computation time seems discouraging. This changes 
if we increase the number of particles with the same mass while keeping the 
same total energy, see Table 2. We also show the dependence of the result 
on the number of collisions which the particles suffer in REGGAE. The 
Nc — 6 case starts to differ slightly from the other algorithms which can 



10 



Table 1 

Mean values of from eq. (17) calculated with five different generators. Number of 
momentum configurations (events) used in the calculations is N. Each event consists 
of n = 5 particles with m = 1 GeV. The last row shows the time for computation 
with 10^ configurations. 



N 


REGGAE (iVc = 6) 


NUPHAZ 


RAMBO 


wGENBOD 


uGENBOD 


10^ 


195.4 


270.48 


181.60 


205.58 


1512 


10^ 


214.5 


223.62 


199.24 


203.81 


1610 


10^ 


211.2 


212.56 


211.63 


208.72 


1635 


10^ 


209.8 


207.71 


208.97 


209.334 


1655 


time 


26 min 


5 min 


1 min 


7 min 


7 min 



Table 2 

Mean values of as in Table 1, but for n = 30 particles with m = I GeV. The 
dependence on the number of collisions in REGGAE is also shown. The last row 
shows the time for computation with 10^ configurations. 



N 


REGGAE 

A^c = 6 


REGGAE 

Nc = 8 


REGGAE 

Nc = 12 


NUPHAZ 


RAMBO 


wGENBOD 


10^ 


13.63 


13.41 


13.78 


12.77 


13.23 


12.23 


10^ 


13.99 


13.52 


13.36 


13.15 


13.21 




10^ 


13.84 


13.42 


13.19 


13.06 


13.12 




time 


16 min 


20 min 


28 min 


6 min 


11 min 


300 min 



Table 3 

Mean values of /s as in Table 1, but for n = 60 particles with m = 1 GeV. The last 
row shows the t ime for computation with 10^ configurations. 



N 


REGGAE {Nc = 12) 


NUPHAZ 


10^ 


0.6339 


0.6185 


10^ 


0.6185 


0.6315 


time for 10^ 


6 min 


63 min 



be fixed by increasing Nc to 12. Our algorithm is clearly most effective in 
cases where the total mass of particles makes up a large part of the energy 
budget in the simulation, see Table 3 where the total mass of particles nm ~ 
60 GeV represents 60% of the energy budget (RAMBO and GENBOD were 
too slow to be included in the Table). This is not surprising: we developped 
REGGAE for the use in simulations of multiparticle production in nuclear and 
high multiplicity hadronic collisions while the other methods were derived for 
applications in high energy physics where most of the available energy often 
goes into momenta. 
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For large n, REGGAE can be the fastest even in the case when 90% of the 
available energy goes into momenta and just 10% is in the particle masses. 
The time to generate 10"^ configurations for n = 500 particles with the mass 
0.1 GeV and the total four-momentum P — (500, 0, 0, 0) GeV/c is 8 min for 
REGGAE {Nc = 10) and 23 min for NUPHAZ. 



5 REGGAE: a reference manual 

The package comes with the following files: 

Mcikef ile for linux-like systems with g++ compiler this is the simple Makefile 
to compile the example program. 

example . cpp is the file with the main routine used to illustrate the action 
of REGGAE. You do not need this file if you want to embed REGGAE in 
your simulation. 

example.dat is the output file of the code if run in the form as it is dis- 
tributed. 

reggae . cpp is the file where main algorithms are implemented. 

reggae . hpp the header file which must be included in the application that 

is supposed to use REGGAE. 

specrel.cpp supporting file with definitions of classes for four-vectors, Lo- 
rentz boosts etc. 

specrel.hpp header file for specrel.cpp. 

To use the method, the user must define the following variables: 
int n the number of particles 

double * amass which must be allocated to include n values for the masses 
of the particles. Before calling the routine the array must be filled with the 

masses. 

vector4 * avec which must be allocated to hold n four-vectors. This array 

will hold the particle four-momenta. 
vector4 P the total four-momentum. Before calling the routines it must 

posses a value. 

long int seed is the initial seed for the random number generator. It may 
be initialised always with the same value or with the help of the machine 
time. 

The four- vector class vector4 is defined in specrel .hpp. It has been designed 
so that its use is intuitive. To get or set the component of the four-vector use 
[*] (e.g. avec[0], numbering runs from to 3). Minkowski metric (-|- ,-,-,-) is 
applied so that a*b gives the proper four- vector product. 
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To generate a sample of momenta in REGGAE, first the GENBOD algorithm 
must be called by 

Mconserv (P , n , amass , avec , feseed) ; 

where the variables have been explained above. Now the momenta are stored 
in avec. The second step is to call 

collision(n,avec,&seed) ; 

This step reshuffles momenta and returns them also in the array avec. This 
array of momenta is the result of the generator. 

The procedure is illustrated in the main routine distributed with the package. 

By default, collisions is set to run eight scatterings per particle. It is pos- 
sible to decrease this number in order to run faster or increase it in order to 

have more confidence that the generated configurations are in the regime of 
saturated entropy. This is done by setting the variable int RG.opai to the 
value equal to half of the number of collisions. 



6 Summciry 

REGGAE gives the possibihty of robust and effective simulation of events 
with few as well as many particles with strict conservation of momentum 
and energy. This is very much needed feature in many models used for the 
simulation of multiparticle production. 
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